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In this paper we propose a tight-binding molecular dynamics with parameters fitted to first- 
principles calculations on the smaller clusters and with an environment correction, to be a powerful 
technique for studying large transition/noble metal clusters. In particular, the structure and stability 
of Cun clusters for n = 3 — 55 are studied by using this technique. The results for small Cu„ clusters 
(n = 3 — 9) show good agreement with ab initio calculations and available experimental results. In 
the size range 10 < n < 55 most of the clusters adopt icosahedral structure which can be derived 
from the 13-atom icosahedron, the polyicosahedral 19-, 23-, and 26-atom clusters and the 55-atom 
icosahedron, by adding or removing atoms. However, a local geometrical change from icosahedral 
to decahedral structure is observed for n = 40 — 44 and return to the icosahedral growth pattern is 
found at n = 45 which continues. Electronic "magic numbers" (n = 2, 8, 20, 34, 40) in this regime 
are correctly reproduced. Due to electron pairing in HOMOs, even-odd alternation is found. A 
sudden loss of even-odd alternation in second difference of cluster binding energy, HOMO-LUMO 
gap energy and ionization potential is observed in the region n ~ 40 due to structural change there. 
Interplay between electronic and geometrical structure is found. 

PACS numbers: 36.40.-c, 36.40.Cg, 36.40.Mr, 36.40.Qv 



I. INTRODUCTION 



The study of clusters has become an increasingly in- 
teresting topic of research in both physics and chemistry 
in recent years, since they span the gap between the mi- 
croscopic and macroscopic materials EllS • Metallic clus- 
ters play a central role in catalysis [ij S HI E| an d nan- 
otechnology 0,0, Ei Cluaters of coinage metals Cu, Ag 
and Au have been used in a wide range of demonstra- 
tion [1 H H H HH 0. The determination of struc- 
tural and electronic properties and the growth pattern of 
coinage metal clusters are of much interest both experi- 
mentally M, El 111 111 BIH E E El and theoret- 
ically |24| . The electronic configura- 
tion of the coinage metals are characterized by a closed d 
shell and a single s valance electron [Cu : Ar(3d) 10 (4 : s) 1 , 
Ag : Kr(Ad) w '(5s) 1 , Au : Ae(5d) 10 (6s) 1 ]. Due to the 
presence of single s electrons in the atomic outer shells, 
the noble metal clusters are expected to exhibit certain 
similarities to the alkali metal clusters. Electronic struc- 
ture of alkali metal clusters are well described by the 
spherical shell model, which has successfully interpreted 
the "magic numbers" in Na„ and K n clusters 0,12 ■ A 
number of experimental features of noble metal clusters 
are also qualitatively well described in terms of simple s 
electron shell model. For instance, the mass abundance 
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spectrum of Cu~, Ag~ and Au~ clusters, which reflects 
the stability of clusters, can be explained by the one- 
electron shell model . But some experimental studies 
[Til IT2I HI Hij indicate that the localized d electrons of 
the noble metals play a significant role for the geomet- 
rical and electronic structure through the hybridization 
with more extended valence s electron . Therefore, it is 
important to include the contribution of 3d electrons and 
the s — d hybridization for Cu n clusters. 

Bare copper clusters in the gas phase have been stud- 
ied experimentally by Taylor et al. fl5| and Ho et al. 
[ill using photoelectron spectroscopy (PES). Knickel- 
bein measured ionization potentials of neutral copper 
clusters and found evidence of electronic shell struc- 
ture [l7^ . Very recently cationic copper clusters have 
been studied using threshold collision-induced dissocia- 
tion (TCID) by Spasov et al. 01- Copper clusters have 
been also investigated theoretically by various accurate 
quantum mechanical and chemical approaches. Masso- 
brio et al. E studied the structures and energetics of 
Cu n (n = 2,3,4,6,8,10) within the local density ap- 
proximation of density functional theory (DF-LDA) by 
using the Car Parinello (CP) method. Calaminici et 
al. [2(j used the linear combination of Gaussian-type 
orbitals density functional (LCGTO-DFT) approach to 
study Cu n , Cu~ and Cit+ clusters with n < 5. Akeby et 
al. [2j| used the configuration interaction (CI) method 
with an effective core potential (ECP) for n < 10. In 
an earlier communication |22l | we studied the small Cu n 
clusters for n < 9 by using full-potential muffin-tin or- 
bitals (FP-LMTO) technique. 
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Ideally, the sophisticated, quantum chemistry based, 
first-principles methods predict both the stable geome- 
tries and energetics to a very high degree of accuracy. 
The practical problem arises from the fact that for ac- 
tual implementation these techniques are limited to small 
clusters only. None of the methods described above 
can be implemented for clusters much larger than ~ 
10 atoms, because of prohibitive computational expense. 
The aim of this communication is to introduce an semi- 
empirical method, which nevertheless retains some of the 
electronic structure features of the problem. The empir- 
ical parameters are determined from first-principles cal- 
culations for small clusters, and corrections introduced 
for local environmental corrections in the larger clusters. 

In recent years empirical tight-binding molecular dy- 
namics (TB-MD) method has been developed as an al- 
ternative to ab initio methods. As compared with ab 
initio methods, the parameterized tight-binding Hamil- 
tonian reduces the computational cost dramatically. The 
main problem with the empirical tight-binding methods 
has always been the lack of transferability of its empiri- 
cal parameters. We shall describe here a technique that 
allows us to fit the parameters of the model from a fully 
ab initio, self-consistent local spin-density approximation 
(LSDA ) ba sed FP-LMTO calculation reported earlier by 
us p3. l23l| for the smaller clusters and then make cor- 
rection for the new environment for clusters in order to 
ensure transferability (at least to a degree). 

It should be mentioned here that Copper clusters 
have also been investigated by other empirical methods. 
D'Agostino carried molecular dynamics using a quasi- 
empirical potential derived from a tight-binding approach 
for nearly 1300 atoms j24|. More recently Darby et al. 
carried geometry optimization by genetic algorithm us- 
ing Gupta potential p5l | for Cu„, Au n and their alloy 
clusters in the size range n < 56 [2(|. These kinds of 
empirical atomistic potentials are found to be good to 
predict ground state geometries but can not predict elec- 
tronic properties such as electronic shell closing effect for 
n =2, 8, 20, 40, highest occupied-lowest unoccupied 
molecular level (HOMO-LUMO) gap energy and ioniza- 
tion potential. Our proposed TBMD scheme will allow 
us to extrapolate to the larger clusters to study both the 
ground state geometries as well as ground state energetics 
as a function of cluster size. 

Menon et al. have proposed a minimal parameter 
tight-binding molecular dynamics (TBMD) scheme for 
semiconductors [27], l29j| and extended the method for 
transition metal (Ni n and Fe n ) clusters [3(3, |^|. Re- 
cently Zhao et al. has applied this method for silver 
clusters . In the present work, we shall introduce a 
similar TB model for copper. 

Using this TBMD method, we shall investigate the 
stable structures, cohesive energies, relative stabilities, 
HOMO-LUMO gaps and ionization potentials of Cu n 
clusters in the size range n < 55. We shall indicate the 
comparison between the present results for small clusters, 
n < 9, with those of our previous FP-LMTO calculations 



and other ab initio and available experimental results. 
This is essential before we go over to the computation- 
ally expensive study of larger clusters. 

II. COMPUTATIONAL METHOD 

Menon et al. introduced a minimal parameter tight- 
binding molecular dynamics (TBMD) scheme for tran- 
sition metal clusters [3(1 [jOJ ■ Here we will describe the 
main ingredients. In this tight-binding scheme the total 
energy E is written as a sum, 

E = E e i + E rep + Ebond- (1) 

E e i is the sum of the one-electron energies for the occu- 
pied states efc, 

occ 

E el = Y,ek, (2) 

fc 

where the energy eigenvalues are calculated by solving 
the eigenvalue equation 

= e*|*fc>, (3) 

where H is the one-electron Hamiltonian and \^k) is 
electronic wave function for fcth level of the eigenstate. 
In the TB formulation, the single particle wavefunc- 
tions \^k) are cast as a linear combination of orthog- 
onalized basis functions 0>i U , in the minimum basis set 

(f = S,P x ,Py,p z ,d x y,dy z ,d zx ,d x 2_y2,d 3z 2_ r 2), 

\*k) = £c* |<M, (4) 

i v 

where i labels the ions. 

The TB Hamiltonian H is constructed within Slater- 
Koster scheme |33| , where the diagonal matrix elements 
are taken to be configuration independent and the off- 
diagonal matrix elements are taken to have Slater-Koster 
type angular dependence with respect to the inter-atomic 
separation vector r and scaled exponentially with the 
inter-atomic separation r: 

V \,W,n = ^A,A' :M ( c 5 '(^ TO ; n ) e ZP[-"( r - rf )L (5) 

where d is the equilibrium bond length for the fee bulk 
copper, S(l,m,n) is the Slater-Koster type function of 
the direction cosines l,m,n of the separation vector r 
and a is an adjustable parameter (= 2/d) |3l| . 

The Hamiltonian parameters are determined from the 
dimensionless universal parameters 77 A A < Q , 

Where is characteristic length for the transition metal 
and the parameter r = for s — s, s—p and p — p interac- 
tions, t = 3/2 for s — d and p — d interactions and r = 3 
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TABLE I: Parameter on site energies, E s , E p and Ed and 
the universal constants 7]\x'u 

for Cu 134 ■ 



TABLE II: The adjustable parameters (f>o, a,b, and c. 



Parameter 


Value 


Parameter 


Value 




0.67 A 


T\pp-E 


-0.81 


E s 


-20.14 eV 


f]sd(T 


-3.16 


E p 


100.00 eV 


Vpda 


-2.95 


Ed 


-20.14 eV 


rjpd-K 


1.36 


f]ssa 


-0.48 


Vdda 


-16.20 


7]spa 


1.84 


Vddn 


8.75 


7]pp(T 


3.24 


TjddS 


0.00 



for d — d interaction. In Table [I] we present the param- 
eter r<2, the on-site energies E s ,E p ,Ed and the universal 
constants rj x x > for Cu [34| ■ According to Ref. [3(J and 
Ref. [U , we set E s = Ed and E p large enough to prevent 
p-orbital mixing |34| . This choice of our tight-binding 
parameters reproduces the band structure of the fee bulk 
Cu crystal given by Harrison |34| . 

The repulsive energy E rep is described by a sum of 
short-ranged repulsive pair potentials, tfiij, which scaled 
exponentially with inter-atomic distance, 

Erep = y ] y ] (f>ij ( r ij ) 
1 j(>i) 

= E E (f> Q exp[-(3(r t:j -d)], (7) 
» j(>«) 

where r»j is the separation between the atom i and j and 
/?(= Aa) is a parameter. E rep contains ion-ion repulsive 
interaction and correction to the double counting of the 
electron-electron repulsion present in E e i. The value of 
0o is fitted to reproduce the correct experimental bond 
length of the Cu dimer 2.22 A ^ is S ivcn in Tabl(jTT| 

The first two terms of the total energy are not sufficient 
to exactly reproduce cohesive energies of dimers through 
bulk structures. Tomahek and Schluter |36j | introduced 
a coordination dependent correction term, Eb on d, to the 
total energy, which does not contribute to the force, it is 
added to the total energy after the relaxation has been 
achieved. However, for the metal clusters, this correction 
term is significant in distinguishing various isomers for a 
given cluster |3lj . 











n 






+ c 











(8) 



where n and nb are the number of atoms and total num- 
ber of bonds of the cluster respectively. Number of bonds 
rib are evaluated by summing over all bonds according to 
cut-off distance r c and bond length 



n b 



E 



exp 



A 



(9) 



The parameters a, b and c in the equation (6) are then 
calculated by fitting the coordination dependent term, 



00 (eV) 


a (eV) 


b (eV) 


c (eV) 


0.034 


-0.0671 


1.2375 


-3.0420 



Ebond, to the ab initio results for three small clusters of 
different sizes according to the following equation 



bond 



E 



ab initio 



Epi — E r 



(10) 



Thus we have four parameters 4>q, a, b, and c in this 
TB model. These parameters are once calculated (given 
in the Table [nj for small clusters to reproduce known 
results (whatever experimental or theoretical) and then 
kept fixed for other arbitrary size cluster. To determine 
the parameters a, b and c we use the experimental bind- 
ing energy of Cu dimer 1.03 eV/atom and the ab 
initio FP-LMTO results for Cu 4 and Cu 6 in Ref. 
For the Cu2 dimer calculated vibrational frequency (226 
cm -1 ) has reasonable agreement with experiment |3^ | 
(265 cm" 1 ). 

In molecular dynamics scheme the trajectories {Rj(t)} 
of the ions are determined by the potential energy sur- 
face E[{Rj(t)}} corresponding to the total energy of the 
electronic system. The force acting on the z-th ion is thus 
given by, 



P 4 = -V Ri E[{Rj}} 



(11) 



This equation can be further simplified by making use 
of the Hellmann-Feynman |38j theorem 



rep • 



(12) 



The second term in the above equation is the short- 
ranged repulsive force. We should note that Pulay cor- 
rection term does not play any role in any semi-empirical 
TBMD [11. The reason is twofold. Within TBMD we 
directly compute the derivative of the TB Hamiltonian 
matrix element and the basis wavefunctions never ap- 
pear explicitely, rather they are implicitely contained in 
the fitted matrix entries. 

The motion of the atoms follow a classical behaviour 
and is governed by the Newton's law : 



d 2 Ki 



(13) 



where m is the atomic mass. 

For numerical simulation of Newtonian dynamics, we 
use the velocity Vcrlct molecular dynamics method 01 
for updating the atomic coordinates, which is given by, 

Ri(t + St) = TLi(t) + Vi(t)6t + -Lf,(<)(^) 2 , (14) 

Am 
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TABLE III: Point group (PG) symmetry, cohesive energy per atom, difference in cohesive energy per atom AE and average 
bond length (r) of the ground state structure and different isomers for Cu n clusters with n < 9 obtained from TB calculation 
and comparison with ab initio calculations fl9l I2H. I2H . AE — 0.00 represents the most stable structure for a particular n. 
Cohesive energy corresponding to the ground state structure in FP-LMTO (2^1. DF-LDA [l9Tl (in parentheses) calculations and 
the values from TCID experiment Hs| are given. For C117, C 3v (I) is the bicapped trigonal biprism and C 3v (II) is the capped 
octahedron. 



Cluster 


PG 

Symmetry 


Present 


Binding Energy (eV/atom) 

Theory" Experiment 6 


Present 


AE (cV/atom) 

Theory c Theory d 


(r) 
(A) 


Cu 3 


C2v 


1.43 


1 60H (S3 1 ) 


1 07+0 1 2 


0.00 


0.00 


2.25 




D 3h 


1.32 






0.11 


0.06 


2.24 






1.13 






0.30 


0.00 


2.24 


Cua 


D 2h 


2.00 


2 00C2 09") 


1.48±0.14 


0.00 


0.00 0.00 


2.23 




D ih 


1.73 






0.27 


0.56 


2.22 




T d 


1.46 






0.54 


0.89 


2.24 


Cu 5 


C2v 


2.24 


2.19 


1.56±0.15 


0.00 


0.00 


2.23 




D 3h 


2.03 






0.21 


0.37 


2.38 


CU6 




2.54 


2.40(2.49) 


1.73±0.18 


0.00 


0.00 0.00 


2.40 






2.40 






0.14 


0.01 


2.39 




o h 


1.98 






0.56 


0.87 0.04 


2.41 


Cu 7 


D 5h 


2.63 


2.65 


1.86±0.22 


0.00 


0.00 


2.41 




C 3v {I) 


2.50 






0.13 


0.32 


2.63 




C 3v (II) 


2.30 






0.33 




2.45 


Cug 


C s 


2.87 


2.73(2.84) 


2.00±0.23 


0.00 


0.20 


2.41 




0,, 


2.64 






0.23 




2.61 






2.57 






0.30 


0.00 


2.59 




T d 


2.51 






0.36 


0.15 


2.39 


CU9 


c> 


2.87 


2.80 




0.00 




2.44 




C2v 


2.84 






0.03 




2.59 




Cs 


2.60 






0.27 




2.41 



"From Kabir et al. fB.ef.l2ll and Massobrio et al. fRef.ligl'). 
b Calculated from Spasov et al. CRef.lls(|'l. 
c From Akeby et al. CRefEHV 
''From Massobrio et al. (Ref.liH'). 



where the velocity Vj of the ith atom at t+8t is calculated 
from Fi at t and t + St as 

Vi(t + St) = Vi(t) + ^-[F t (t) + Fi(t + St)]St. (15) 
2m 

At this stage most authors carry out either dissipative 
dynamics or free dynamics with feedback ^lj. The rea- 
son for this is as follows : for numerical integration of 
Newton's equations we have to choose a finite time-step 
St. Ideally this should be as small as possible, but that 
would require an excessively long time for locating the 
global minimum. However, a large choice of St leads to 
unphysical heating up of the system, leading to instabil- 
ity. Dissipative dynamics has been suggested as a way 
of overcoming this. We add a small extra friction term 
carefully F F — 7mR |3lJ ■ In the present calculation 
7m = 0.32 amu/psec, and the time step St is taken to 
be 1 fsec and the total time for molecular dynamics sim- 
ulation is ~ 100 - 200 pscc, depending upon the cluster 
size and initial cluster configuration with the several an- 
nealing schedule. Methfessel and Schilfgaarde ^l| have 
also used an alternative technique of free dynamics with 
feedback to overcome the above difficulty. 

The results of the molecular dynamics may depend sen- 
sitively on the starting configuration chosen. The final 



equilibrium configurations often correspond to local min- 
ima of the total energy surface and are metastable states. 
For the smaller clusters simulated annealing can lead to 
the global minimum. We have found the global mini- 
mum configurations of the smaller clusters by the simu- 
lated annealing technique. However, this is often not the 
case for the larger clusters. Recently more sophisticated 
techniques like the genetic algorithm has been proposed 
■ We have not tried this out in this communica- 
tion, but propose this as an efficient technique for further 
work. 



III. RESULTS AND DISCUSSION 

A. Geometry optimization 

We have applied this TBMD scheme to Cu n clusters 
for n < 55. Since the present scheme imposes no a priori 
symmetry restrictions, we can perform full optimization 
of cluster geometries. For small clusters (n < 9) we can 
able to perform a full configurational space search to de- 
termine the lowest-energy configuration. Here they serve 
as a test case for the calculation of larger clusters with 
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FIG. 1: Most stable structures for copper clusters with n — 10 — 55 atoms. Most of the clusters adopt icosahedral structures 
except for n = 40 — 44, where the structures are decahedral. 



n > 10. In Table ITTT1 we present a detailed comparison 
of binding energy per atom, difference in binding energy 
AE and average bond length (r) for n < 9 with available 
experimental [3 and ab initio 0, l2lTl2^| results. We 



found, i n agree ment with experimental |16( and theoret- 
ical [Ial20il23 results, very small copper clusters (Cu 3 , 
Ciii and C115) prefer planer structures. More detailed 
comparison, with experimental and ab initio results, can 
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be found elsewhere |46[ . 

From the present results and detailed comparisons with 
various experimental 0,0] and ab initio I20I I2H I22I 
E3jEME3 results available, we find reasonable agreement 
among this TBMD scheme and ab initio calculations for 
small clusters with n < 9 0], which allow us to con- 
tinue the use of this TBMD scheme for larger clusters 
with n > 10. For larger clusters (10 < n < 55), due 
to increasing degrees of freedom with cluster size, a full 
configurational search is not possible with the available 
computational resources. Instead, led by the experimen- 
tal and theoretical results on small clusters, we examined 
structures of various symmetries for each size. Most sta- 
ble structures for n = 10 — 55 atom clusters are given in 
FigU] 

In this regime, the structures predicted by this TB 
model are mainly based on icosahedron. The most sta- 
ble structure of Cuy is a pentagonal bipyramid (D§h 
symmetry; see Tabl dnT|) . which is the building block for 
the larger clusters with n > 10. For Ouio, wc found a 
tricapped pentagonal bipyramid to be the most stable 
structure. Ground state structures of Cuu and Cuyi 
are the uncompleted icosahedron with lack of one and 
two atoms respectively and a Jahn- Teller distorted first 
closed shell icosahedron is formed at Cuv$. For Cui 3 , the 
fee like cuboctahedron is less stable than the icosahedron 
by an energy 0.05 eV per atom. In agreement Lammcrs 
and Borstel, on the basis of tight-binding linear muffin- 
tin orbital calculations, was also found the icosahedron 
to be the ground state of Cu\s, though the difference 
in energy between the icosahedron and the cuboctahe- 
dron was calculated to be only 0.2 eV/atom 0- The 
ground state structures for Cuu, Cuis, Cu\q, and Cun 
are the 13- atom icosahedron plus one, two, three and 
four atoms respectively. A double icosahedron is formed 
for Cwig {D$h symmetry). This structure has two inter- 
nal atoms, 12 six-coordinate atoms at either end and five 
eight-coordinate atoms around the waist of the cluster. 
Based on the structure for Cu\g, the stable Cuis cluster 
is a double icosahedron minus one of the six-coordinate 
atoms located at either end (Csu symmetry). Icosahe- 
dral growth continues for 20 < n < 55 atom clusters. 
Polyicosahcdral structure in the form of a " triple icosa- 
hedron" (Dsh symmetry; the structure can be viewed 
as three interpenetrating double icosahedra) is the most 
stable structure for C1123 cluster. The next closed shell 
polyicosahedra is found for Cu2e cluster. Finally, the 
second closed shell icosahedron is formed for CU55 which 
is more stable than the closed cuboctahedral structure 
by an energy difference 6.27 eV. This can be explained 
in terms of their surface energy. The surface energy 
of the icosahedral structure is lower than that of the 
cuboctahedral structure, because the atoms on the sur- 
face of the icosahedron are five-coordinate compared to 
the four-coordinate atoms on the surface of the cuboc- 
tahedron. In our calculation, exception to the icosahe- 
dral growth is found at around Cu^o- The situation re- 
garding geometrical structure in this size range is more 



complex. The structures for n = 40 — 44 atom clus- 
ters are oblate, decahedron like geometries. Return to 
the icosahedral structure is found at n = 45. In the size 
range n = 40 — 44, the structural sequence is decahedron- 
icosahedron-cuboctahedron in decreasing order of stabil- 
ity, whereas in the region n = 45 — 55, the structures 
retain icosahedron-decahedron-cuboctahedron sequence. 

This results are in agreement with the experimental 
study of Winter and co-workers [5(j , where they found 
a bare copper cluster mass spectrum recorded with ArF 
laser ionization shows a sudden decrease in the ion sig- 
nal at Cu~±2 > an d from this observation they argued that 
a change in geometrical structure might occur there, 
though they have not concluded about the nature of this 
geometrical change. They also found a dramatic decrease 
in water binding energy for Cu^q and Cu^i, and con- 
cluded that this may represent a return to the icosahedral 
structure as the closed shell at CU55 is approached. 

D'Agostino [24[ confirms our prediction, who per- 
formed molecular dynamics using a tight-binding many- 
body potential and found that icosahedral structures are 
prevalent for clusters containing less than about 1500 
atoms. Valkealahti and Manninen [Hi]], using effective 
medium theory, also found icosahedral structures are en- 
ergetically more favourable than the cuboctahedral struc- 
tures for sizes up to n ~ 2500 is consistent with our re- 
sult: Fig|3| shows cuboctahedral structures are least sta- 
ble among the three structures, icosahedron, decahedron 
and cuboctahedron. By contrast, Christcnscn and Ja- 
cobsen jF^ predicted more fcc-like structures in the size 
range n = 3 — 29, in their Monte Carlo simulation us- 
ing an effective medium potential. But they corr ectly 
reproduced the "magic numbers" in that regime [^3, l53j |. 

These results can be compared with the genetic algo- 
rithm study on copper clusters by Darby et al. (2(| , using 
Gupta potential. In agreement with the present study, 
Darby et al. found that most of the clusters in this regime 
adopt structures based on icosahedron. They also found 
exceptions to the icosahedral growth at around Cu^o, 
where the structures adopt decahedron li ke g eometries 
(exact numbers are not available in the Ref.[26|). But the 
present study disagree with the genetic algorithm study 
in two points. Firstly, for 25 atom cluster, they found a 
more disordered structure, while the present study pre- 
dict it to be an icosahedron based structure which can 
be derived by removing one surface atom from the 26- 
atom polyicosahcdron. Finally, they found an fcc-like 
truncated octahedral structure for Cu^. Instead, the 
present study predict the icosahedron based structure to 
be the ground state, where this structure is energetically 
more favourable than the truncated octahedral structure 
by an energy AE = 0.17 eV/atom. Although the ge- 
netic algorithm search for global minima is more efficient 
technique than molecular dynamics, use of the empirical 
atomistic potential is the main reason |54j] for this kind 
of disagreement between Darby et al. and the present 
study. 



B. Binding energies and relative stabilities 

The computed size dependence of the binding energy 
per atom for Cu n clusters with n = 2 — 55 is depicted 
in Fig|21 (upper panel). Among all the isomeric geome- 
tries examined for a certain cluster size n, the highest 
cohesive energy has been considered for the FigEI The 
overall shape of the curve matches the anticipated trend: 
binding energy grows monotonically with increasing the 
cluster size. Inset of the Fig|3 (upper panel) shows the 
comparison of our calculated binding energy with the ab 
initio and experimental 0] results. Experimen- 

tally the binding energies of the neutral clusters were 
derived from the dissociation energy data of anionic clus- 
ters from the TCID experiment IhSf and using electron 
affinities from the PES experiment [l6j- The inset figure 
shows that our calculated binding e nerg ies are in good 
agreement with those from DF-LDA |l9| and our previ- 
ous FP-LMTO calculations. However, our binding 
energies are systematically overestimated, by an energy 
0.53 ± 0.12 to 0.79 ± 0.22, than the experimental bind- 
ing energies. The LDA based ab initio calculations al- 
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FIG. 2: (Upper panel) Binding energy per atom as a function 
of cluster size n 1//3 . Inset of the upper panel represents a com- 
parison of binding energy per atom as a function of cluster size 
n, among the present TBMD (□), FP-LMTO (O), DF-LDA 
(A) calculations and experimental (o) values. (Lower panel) 
Variation of relative stability A2E with cluster size n. Shell 
closing effect at n = 8, 18, 20, 34, 40 and even-odd alternation 
up to n ~ 40 are found. However, due to geometrical effect 
this even-odd alternation is disturbed at n = 11, 13 and 15. 
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FIG. 3: Comparison of binding energies per atom as a func- 
tion of cluster size n among cuboctahedral, decahedral and 
icosahedral structures. For the whole region most of the clus- 
ters prefer icosahedral structure. However, a local geometrical 
change from icosahedral to decahedral structure is found for 
n = 40 - 44. 



ways over-estimate binding energies. This is a character- 
istic of the LDA. In the present study, TB parameters 
have been fitted to the ab initio LDA calculations for 
very small calculations ■ It is not surprising therefore 
that the binding energies are over-estimated. In fact, 
the present results agree well with other LDA based cal- 
culations 0, , all of which overestimate the binding 
energy. 

In the FigOll we compared binding energy per atom 
for cuboctahcdral, decahedral and icosahedral structures 
for the clusters containing n = 30 — 55 atoms. Fig|3| 
shows most the clusters in this size range have icosa- 
hedral structures. However, a local structural change 
occured for n = 40 — 44, where the structures adopt dec- 
ahedral structure rather than icosahedral one. Return 
to the icosahedral growth pattern is found at n = 45 
and continues up to 55-atom cluster. From the Fig[3]it is 
clear that among cuboctahedral, decahedral and icosahe- 
dral structures, cuboctahedral structures are least stable 
than the other two. 

The second difference in the binding energy may be 
calculated as 



A 2 E(n) = E{n + 1) + E{n - 1) - 2E(n), 



(16) 



where E(n) represents the total energy for an n-atom 
cluster. A2E(n) represents the relative stability of an n- 
atom cluster with respect to its neighbors and can be di- 
rectly compared to the experimental relative abundance 
: the peaks in AzEin) coincide with the discontinuities 
in the mass spectra. These are plotted in the lower panel 
of FigEl We found three major characteristics in the 
Fig|21 (lower panel ). Firstly, even-odd (even > odd) 
oscillation is found. This can be explained in terms of 
electron pairing in HOMOs. Even (odd) clusters have an 
even (odd) number of electrons and the HOMO is dou- 
bly (singly) occupied. The electron in a doubly occupied 
HOMO will feel a stronger effective core potential be- 
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cause the electron screening is weaker for the electrons 
in the same orbital than for inner shell electrons. Thus 
the binding energy of the valence electron with an even 
cluster is larger than of an odd one. This even-odd al- 
ternation is prominent up to n ~ 40. Secondly, due to 
electronic shell or subshell closing, we found particular 
high peak for n = 8, 18, 20, 34 and 40. Unfortunately, 
the present study does not show any evidence of shell 
closing for Cu2 in A2E(n). Finally, the even-odd al- 
ternation is reversed for n = 10—16 with maxima at 
Cuii, CU13 and Cuig, which manifests the geometrical 
effect and therefore there is no peak at n = 14 due to 
electronic subshell closing. Simultaneous appearance of 
these three features in A2-E/(n) demonstrates the inter- 
play between electronic and geometrical structure, which 
is in a gree ment with the experimental study of Winter 
et al. [50|. They found both jellium-like electronic be- 
haviour and icosahedral structure in copper clusters. In 
an experimental study of mass spectra of ionized copper 
clusters substantial discontinuities in mass spectra 
at n = 3, 9, 21, 35, 41 for cationic and n = 7, 19, 33, 39 for 
anionic clusters as well as dramatic even-odd alternation 
are found. From the sudden loss in the even-odd alterna- 
tion at CM42 in the KrCl spectrum, Winter et al. argued 
about the possible geometrical change there. Therefore, 
we conclude in the section that sudden loss in the A2-E vs 
n plot (lower panel of Fig|2J) is due the structural change 
in that regime. 

Such kind of electronic effects can not be reproduced 
by empirical atomistic potentials. Darby et al. j2||, us- 
ing the Gupta potential, found significant peaks at n = 
7, 13, 19, 23 and 55 due to icosahedral (or polyicosahe- 
dral) geometries. In the present study, we have found 
a peak at n = 13, but not at the other sizes found by 
them. However, the stable structures predicted by us 
are the same : the lowest energy structure of Cu-j is a 
pentagonal bipyramid (D^h symmetry) ; for n = 13 and 
55, the structure are the first and second closed icosahe- 
dral geometries respectively. Polyicosahcdral structures 
are found for n = 19 (double icosahedron) and n = 23 
(triple icosahedron) atom clusters. As the reason, the 
present study shows significant high peaks at Cug, Cui 8 
and Cw,2o due to electronic shell closing effect and av- 
erage peaks at CU22 and Cu^ due to electron pairing 
effect. At these sizes, the electronic effects dominates 
over the geometrical effects and consequently the above 
peaks cannot observed by Darby et. al.. 



C. HOMO-LUMO gap energies 

Besides the second difference of the cluster binding 
energy, a sensitive quantity to probe the stability is 
the highest occupied-lowest unoccupied molecular level 
(HOMO-LUMO) gap energy. In the case of magic clus- 
ters shell or subshell closer manifests themselves in par- 
ticularly large HOMO-LUMO gap, which was previ- 
ously demonstrated experimentally [l(| |55j . Calculated 
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FIG. 4: Highest occupied - lowest unoccupied molecular or- 
bital (HOMO-LUMO) gap energy vs cluster size n. Electronic 
shell closer at n — 2, 8, 18, 20, 34, 40 and even-odd alternation 
are observed. However, sudden loss in even-odd alternation 
is found around n ~ 40 due to the structural change there. 

HOMO-LUMO gap energies are plotted in the FigH 
where we observed even-odd alternation due to electron 
pairing effect and particularly large gap for n = 2, 8, 
18, 20, 34 and 40 due to electronic shell and subshell 
closing. However, sudden loss of even-odd alternation is 
found around n ~ 40 due to the change in the geometrical 
structure in that region. Winter et al. [5(j also found a 
sudden loss in even-odd alternation in the KrCl spectrum 
at Cu42 and concluded this may coincide with any pos- 
sible change in t he g eometrical structure there. In fact, 
Katakuse et al. [lfj observed identical behaviour in the 
mass spectra of sputtered copped and silver cluster ions : 
a dramatic loss of even-odd alternation at n = 42, signi- 
fying a sudden change to a geometrical structure in which 
stability, and abundance, is less sensitive to electron pair- 
ing. Therefore, the sudden loss in the FigQ] again con- 
firms the structural change there. So, the present study 
correctly predicts the "magic numbers" in this regime 
correctly and confirms the experimental prediction : a 
geometrical change (icosahedron — > decahedron) is oc- 
curring around n ~ 40. 



D. Ionization potentials 

Within the present TB scheme, we can get a 'quali- 
tative' description of the ionization potentials with clus- 
ter size according to Koopmans' theorem. This limita- 
tion arises mainly from the choice of the Slater-Kostcr 
(SK) TB parameters and the extent of their transferabil- 
ity [561 ] , which may be improved by the proposed scaling 
scheme of Cohen, Mehl and Papaconstantopoulos • 
However, our aim is to get only a qualitative description 
of ionization potential with cluster size. Calculated ion- 
ization potentials are plotted in the Fig|5] In fact, we 
observed same pattern as it was in HOMO-LUMO gap 
energy vs cluster size : peaks at n = 2, 8, 18, 20, 34, 
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Cluster Size n 

FIG. 5: Ionization potential vs cluster size n. Electronic 
shell closing effect and prominant even-odd alternation up to 
n ~ 40 are observed. 

40 and even-odd alternation due to the same reasons dis- 
cussed in the Sec. IIII Bl and Sec. IIII CI Sudden loss in 
even-odd alternation around n ~ 40 is again confirmed 
from the Fig0 which is due to the geometrical change 
there. 



IV. CONCLUSION 

Using tight-binding model we calculated ground state 
geometries, binding energies, second differences in bind- 
ing energy, HOMO-LUMO gap energies and ionization 
potentials for copper clusters in the size range 2 < 
n < 55. We have fitted the parameters of the present 
TB scheme from our previous ab inito calculations |22j |. 
For small clusters n < 9, present results show good 
agreement with exp erimental [l^ . Il8j | and theoretical 
US, HI El El M, 13 E3 results, which allow us to 
go over the larger size range, 10 < n < 55. 

In the size range 10 < n < 55 most of the clus- 
ters adopt icosahedral geometry which can be derived 
from the 13- atom icosahedron, the polyicosahcdral 19-, 
23-, and 26-atom clusters and 55-atom icosahedron, by 
adding or removing atoms. However, exceptions to the 
icosahedral growth is found around n ~ 40. A local geo- 
metrical transition is found for n = 40 — 44- atom clus- 



ters. This is in agreement with the prediction of the two 
experimental studies by Katakuse et aL^(j and Winter 
et aZ.|H3) where they predicted that a local geometrical 
transition may occur at n = 42, though their results are 
not decisive about the nature of this geometrical change. 
Present results show that around n ~ 40 structures 
are changing from icosahedral to decahedral structure, 
where the structural sequence is dccahcdron-icosahedron- 
cuboctahedron in the decreasing order of stability. Re- 
turn to the icosahedral growth is found at n = 45, with 
the sequence icosahedron-decahedron-cuboctahedron in 
the decreasing order of stability. 

As we have fitted the parameters of the present TBMD 
scheme from LDA based ab inito calculations , calcu- 
lated binding energies are in good agreement with the 
LDA based ab initio calculations but overestimates the 
same calculated from the TCID experiment [l^. In the 
present scheme, the "magic nimbcrs" (n = 2, 8, 18, 20, 
34 and 40) due to electronic shell and subshell closing 
are correctly reproduced in the studied regime. Second 
difference of binding energy, HOMO-LUMO gap energy 
and ionization potential show even-odd oscillatory be- 
haviour because of electron pairing in HOMOs in agree- 
ment with experiment. However, a sudden loss in even- 
odd alternation is found around n ~ 40 in the variation of 
second difference in binding energy, HOMO-LUMO gap 
energy and ionization potential with cluster size. This is 
in agreement with the experimental studies 0, |5(| . We 
conclude this is due to the geometrical change (icosahe- 
dron — > decahedron) around there. Present results show 
that electronic structure can coexist with a fixed atomic 
packing. 

Due to lower computational expense this TBMD 
scheme, with parameters fitted to first-principle calcu- 
lation for the smaller clusters and with an environment 
correction, is a very efficient technique to study larger 
clusters, particularly with n > 10. 
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